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EXECUTIVE  SUMMARY 


This  is  the  final  report  for  the  SERDP  project  MM-1640  and  it  covers  the  research 
results  accomplished  since  the  inception  of  the  project  in  2007.  The  basic  premise  of 
this  project  is  the  theoretical  understanding  and  algorithm  development  of  principal 
component  analysis  (PCA)  as  a  de-noising  and  signal-separation  tool  for  transient 
electromagnetic  (TEM)  data  in  unexploded  ordnance  (UXO)  applications.  There  is 
an  express  need  for  techniques  to  reduce  the  presence  of  random  noise  in  TEM  data 
as  well  as  reduce  the  influence  of  correlated  noise  clue  to  a  wide  variety  of  sources 
on  automatic  anomaly-picking  routines  for  more  accurate  detection  with  fewer  false 
anomalies.  We  have  developed  an  algorithm  and  workflow  for  the  processing  and 
inversion  of  TEM  data  that  attenuates  signal  from  undesired  sources  and  accurately 
inverts  TEM  data  for  diagnostic  UXO  parameters. 

The  research  performed  over  the  life-span  of  this  project  has  progressed  satis¬ 
factorily  with  the  major  research  tasks  completed  approximately  to  project  plan.  In 
addition,  we  have  identified  an  additional  area  of  research  critical  to  the  application 
of  PCA  to  TEM  data  and  added  these  aspects  to  the  research  plan. 

First,  we  developed  a  principal  component  analysis  algorithm  tailored  to  un¬ 
exploded  ordnance  applications.  Decay  characteristics  of  TEM  data  preclude  the 
standard  Karhunen-Loeve  transform;  we  have  addressed  these  issues  with  algorithm 
modifications  and  incorporated  these  into  the  workflow. 

Secondly,  we  identified  the  optimum  choice  of  principal  components  for  the  at¬ 
tenuation  of  both  random  noise  and  correlated  noise,  leaving  the  signal  due  to  UXO 
intact.  We  show  that  the  processed  data  is  optimally  prepared  for  automatic  anomaly 
picking  routines  with  a  highly  reduced  number  of  false  anomalies.  We  demonstrate 
this  on  both  synthetic  examples  of  UXO  surveys,  as  well  as  on  TEM  data  from 
Kaho’olawe,  Hawaii,  USA. 

iii 


Finally,  we  have  identified  a  critical  issue  with  inversion  of  processed  data  that 
results  in  extremely  inaccurate  recovered  models  without  the  incorporation  of  the 
PCA  process  into  the  forward  model.  We  developed  an  inversion  algorithm  which 
takes  the  processing  steps  into  account  during  construction  of  the  inverse  kernels. 
This  leads  to  more  accurate  recovered  models  of  inverted  anomalies. 
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CHAPTER  1 


INTRODUCTION 


1.1  SERDP  relevance 

This  project  addresses  the  need  for  improved  signal  processing  outlined  in  the 
Statement  of  Needs  MMSEED-08-01,  and  has  developed  new  technologies  for  enhanc¬ 
ing  time-domain  electromagnetic  data  (TEM)  by  estimating  and  removing  undesirable 
components  such  as  random  noise  and  other  responses  unrelated  to  metallic  objects. 
The  goal  is  to  increase  the  signal-to-noise  ratio  of  the  TEM  data  as  well  as  to  gen¬ 
erate  an  estimate  of  data  noise  characteristics  (statistical  distribution  and  associated 
parameters)  for  use  in  subsequent  inversion-based  discrimination. 

One  of  the  most  effective  geophysical  techniques  for  UXO  application  is  the  tran¬ 
sient  electromagnetic  method.  Using  a  loop  as  a  source,  a  time-varying  magnetic  field 
is  generated  at  the  surface  which  in  turn  induces  electrical  currents  in  the  ground,  and 
more  importantly  in  buried  metallic  objects.  Once  induced,  these  currents  dissipate 
over  time  due  to  ohmic  losses,  leading  to  a  measured  transient  decay  in  magnetic 
field  flnx  density  at  the  surface  by  geophysical  sensors.  The  currents  induced  in  the 
surrounding  medium  and  in  the  buried  metallic  objects  decay  at  different  rates  and 
can  be  separated.  For  metallic  objects,  the  decay  rates  are  directly  related  to  the 
size,  shape,  depth,  and  electrical  conductivity  of  the  target  body.  Unfortunately,  real 
TEM  data  for  UXO  detection  and  discrimination  are  contaminated  with  outside  noise 
sources.  These  sources  of  noise  can  be  spatially  correlated,  such  as  near-surface  geol¬ 
ogy,  micro-topography,  coil  orientation,  sensor  motion,  and  positioning  errors;  or  they 
can  be  uncorrelated,  such  as  random  RF  interference,  telluric  sources,  and  instrument 
interference.  Each  of  these  sources  of  noise,  in  addition  to  buried  UXO,  contributes 
to  the  total  response  measured  by  the  sensors.  The  effect  of  noise  is  especially  strong 
at  late-times  in  the  TEM  data.  Noise  from  such  numerous  and  varying  sources  often 
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limits  the  ability  of  the  TEM  method  to  adequately  discriminate  hazardous  muni¬ 
tions  from  non-hazardous  items.  Thus,  there  is  a  need  for  developing  methods  to 
enhance  TEM  data  in  UXO  discrimination,  and  the  key  to  the  solution  under  such 
conditions  is  to  develop  a  practical  algorithm  that  can  reliably  identify  and  remove 
these  individual  sources  of  noise  component  by  component. 

Of  the  two  types  of  noise,  correlated  noise  is  more  difficult  to  isolate  and  remove 
from  the  data  than  uncorrelated  noise.  Current  methods  of  noise  analysis  in  UXO 
surveys  include  simply  thresholding  a  noise  level  and  ignoring  any  signal  below  the 
chosen  value  (Pasion  and  Oldenburg,  2001b),  stacking,  and  median  filters  to  de-trend 
the  data  (Pasion  and  Oldenburg,  2001a).  To  date,  there  has  not  been  a  concentrated 
research  effort  focused  on  separating  the  various  sources  of  noise,  and  identifying  the 
effects  that  these  individual  sources  have  on  the  data. 

Principal  component  analysis,  or  PCA,  is  a  technique  for  simplifying  a  data 
set  by  reducing  multidimensional  data  sets  to  lower  dimensions  for  analysis.  PCA 
is  basically  an  orthogonal  linear  transformation  that  transforms  the  data  to  a  new 
coordinate  system  such  that  the  greatest  variance  by  any  projection  of  the  data 
comes  to  lie  on  the  first  coordinate  (called  the  first  principal  component),  the  second 
greatest  variance  on  the  second  coordinate,  and  so  on.  By  properly  analyzing  the 
various  components,  one  can  isolate  the  individual  correlated  signals  not  associated 
with  the  UXO,  and  remove  those  components  from  the  TEM  data.  In  an  ideal  case, 
this  leaves  only  the  signal  from  UXO  and  UXO-like  anomalies,  eliminating  completely 
any  error  in  the  signal  caused  by  external  noise. 

We  have  developed  a  set  of  processing  techniques  based  on  the  PCA  algorithm  for 
TEM  data  in  UXO  applications  as  well  as  the  parametric  inversion  of  such  processed 
data  for  recovering  both  extrinsic  and  intrinsic  parameters  for  use  in  discrimination. 
These  techniques  can  be  combined  with  existing  workflows  to  rapidly  identify  poten¬ 
tial  UXO  targets  for  further  processing  with  minimal  required  user  input. 
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1.2  Technical  objectives 


The  purpose  of  this  study  was  to  develop  a  practical  algorithm  that  enhances  the 
signal  in  transient  electromagnetic  data  for  UXO  applications.  We  have  developed 
a  de-noising  algorithm  that  not  only  removes  random,  uncorrelated  noise,  but  also 
separates  the  signals  dne  to  various  sources  through  principal  component  analysis. 
We  separated  the  research  into  aspects  of  PCA  that  address  uncorrelated  signals  and 
aspects  that  address  correlated  but  undesired  signals.  Specifically,  this  study  had  six 
objectives  in  three  categories: 

•  Method  development  for  uncorrelated  noise 

—  Develop  a  PCA  algorithm  tailored  to  TEM  data  from  UXO  surveys,  and 

—  Produce  a  stable  and  automated  algorithm  for  removal  of  uncorrelated 
noise. 

•  Method  development  for  correlated  noise 

—  Study  the  physical  connection  between  signals  due  to  different  sources  and 
the  components  recovered  from  PCA 

—  Determine  the  feasibility  of  using  PCA  to  reduce  TEM  data  in  UXO  appli¬ 
cations  to  the  signal  exclusively  produced  by  UXO  and  UXO-likc  anoma¬ 
lies,  and 

—  Apply  the  developed  algorithm  to  data  examples  and  assess  the  effective¬ 
ness  through  discrimination/classification  tests. 

•  Method  development  for  inversion  of  processed  data 

—  Develop  a  method  for  parametric  inversion  of  processed  data  for  recovery 
of  diagnostic  parameters  in  UXO  applications. 
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1.3  Project  summary 


In  this  report,  we  summarize  the  major  research  accomplishments  of  the  SERDP 
MM-1640  project.  In  short,  we  have  accomplished  the  main  objectives  of  this  project, 
leading  to  a  workflow  ready  to  be  inserted  into  industry-standard  practices  for  the 
improved  detection  of  UXO  in  noisy  environments. 

In  Chapter  2  we  begin  by  presenting  background  material  on  principal  compo¬ 
nent  analysis  which  will  be  used  throughout  this  report.  We  first  develop  necessary 
modifications  to  the  basic  PCA  routines  in  Chapter  3  before  presenting  de-noising 
algorithms  for  both  uncorrelated  and  correlated  sources  of  noise.  We  continue  in 
Chapter  3  by  discussing  requirements  and  procedures  for  large-scale  PCA  decompo¬ 
sitions.  We  conclude  the  chapter  by  showing  difficulties  in  inverting  PCA-processed 
data  and  develop  an  algorithm  for  both  parametric  and  non-parametric  inversion  of 
this  processed  data. 

In  Chapter  4,  we  provide  a  series  of  examples  demonstrating  the  efficacy  of  PCA 
in  transient  electromagnetic  survey  processing.  We  begin  by  showing  how  PCA  can 
remove  uncorrelated  noise  in  a  field  dataset  from  Kaho’olawe,  HI.  We  then  present  a 
parametric  inversion  result  showing  that  PCA  can  improve  the  recovery  of  diagnostic 
parameters  for  UXO/scrap  metal  discrimination.  We  conclude  the  examples  by  taking 
a  dataset  contaminated  by  correlated  noise  through  a  complete  workflow  in  UXOLab 
both  with  and  without  processing  by  PCA.  We  show  a  marked  reduction  in  the 
number  of  false  anomalies  without  compromising  the  detection  of  actual  UXO. 

We  conclude  this  report  in  Chapters  5  and  6  with  a  discussion  of  the  project 
results,  recommendations,  and  publications/presentations  generated  throughout  the 
lifespan  of  SERDP  Project  MM-1640. 
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CHAPTER  2 


BACKGROUND 


2.1  Principal  component  analysis 

Principal  component  analysis  (PCA)  is  a  statistical  method  for  analyzing  ob¬ 
servations  in  multiple  channels  through  an  orthonormal  projection.  Not  only  does 
principal  component  analysis  have  the  ability  to  separate  and  remove  uncorrelated 
noise,  but  also  to  decompose  a  signal  into  constituent  components  from  its  sources 
in  many  cases.  Consequently,  PCA  has  been  used  for  many  different  applications 
including  digital  image  enhancement,  facial  recognition,  data  transmission  and  com¬ 
pression  (Jones  and  Levy,  1987),  de-noising  radiometric  data  (Minty  and  Hovgaard, 
2002),  and  seismic  de-noising  (Jones  and  Levy,  1987;  Jackson  et  ah,  1991),  to  name 
a  few.  Recently,  PCA  has  been  applied  to  airborne  EM  (AEM)  surveys  to  construct 
intuitive  RGB  maps  based  on  the  principal  components  (Green,  1998).  In  other 
electromagnetic  applications,  Asten  (2009)  has  used  PCA  in  unexploded  ordnance 
TEM  surveys  to  classify  targets,  while  Hu  and  Collins  (2004),  Hu  et  ah  (2004),  and 
Throckmorton  et  ah  (2007)  have  used  similar  techniques  using  higher-order  moments 
(Independent  Component  Analysis)  to  recover  diagnostic  parameters  of  differing  ord¬ 
nance  items.  The  most  common  use  of  Principal  Component  Analysis  in  EM  is  to 
combine  magnetic  and  electromagnetic  datasets  (e.g.  Rose-Pehrsson  et  ah,  1998). 

Although  PCA  methods  vary  greatly  in  their  application,  all  linear  PCA  al¬ 
gorithms  are  similar  in  that  they  deconstruct  a  multi-channel  signal  into  a  set  of 
orthonormal  bases  of  decreasing  energy.  These  sets  can  be  reconstructed  into  the 
original  signal  exactly,  or  a  truncated  set  can  be  used  in  reconstruction  that  tends 
to  eliminate  noise.  Correlated  signals  are  typically  reconstructed  by  added  the  bases 
(from  highest  to  lowest  energy)  until  the  desired  energy  level  is  reached-generally 
chosen  through  a  statistical  analysis  of  error. 
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In  most  surveys  with  multiple  stations  that  record  traces  (decay  curves)  of  data, 
there  are  two  organizations  over  which  a  signal  can  be  said  to  be  “correlated.”  One 
can  analyze  the  correlation  between  time-slices  or  between  traces.  In  essence,  we 
are  either  analyzing  similarities  in  anomaly  geometry  or  similarities  in  decay.  Which 
correlation  we  analyze  will  have  important  implications  on  computation  cost  of  the 
calculations. 

Geometrically,  the  data  set  from  a  TEM  survey  can  be  envisioned  as  a  cloud  of 
'm  x  n‘  measurements  in  an  n-dimensional  data  space  (where  m  is  the  number  of 
channels  and  n  the  number  of  stations,  for  purposes  of  this  derivation).  The  data 
coordinate  system  is  rotated  until  the  first  axis  is  along  the  direction  where  the  data 
have  the  greatest  variance.  The  second  coordinate  axis  is  chosen  to  have  the  next 
highest  data  variance  subject  to  the  constraint  that  it  must  be  orthogonal  to  the 
first,  and  so  on  (Green,  1998).  This  is  accomplished  by  using  the  eigenvectors  of  the 
covariance  matrix  as  the  new  bases,  as  they  will  align  with  these  directions  of  variance 
and  will  be  mutually  orthogonal. 

To  define  these  principal  components,  we  introduce  the  Karhunen-Loeve  Trans¬ 
form  as  a  linear  PGA  tool  as  in  Jones  and  Levy  (1987).  However,  any  eigenvector 
or  singular  value  decomposition  will  produce  an  equivalent  result  to  the  accuracy  of 
that  decomposition  (Minty  and  McFadden,  1998). 

Assume  a  set  of  data  consisting  of  multiple  traces  each  with  multiple  data  values 
Xi(tj),  i  =  1, ...,  n;  j  =  1,  ...m 

where  the  i’th  trace  is  associated  with  a  particular  location  and  the  index  j  corre¬ 
sponds  to  a  particular  channel  of  observation  at  that  given  location.  For  example, 
one  may  have  an  EM63  instrument  occupying  n  locations  along  a  line  or  within  a 
grid  during  a  UXO  survey,  which  would  yield  n  traces  of  decaying  voltage,  each  con¬ 
taining  25  channels  of  data  (m  =  25).  As  stated  earlier,  principal  component  analysis 
first  generates  a  rotation  matrix  that  rotates  the  data  traces  onto  a  set  of  orthogonal 
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directions  called  principal  component  directions.  In  such  a  rotation,  the  resulting 
components  are  ordered  by  decreasing  energy.  The  earlier  components  with  most 
energy  tend  to  capture  the  coherent  data  signal  whereas  the  later  components  tend 
to  represent  incoherent  noise  in  the  entire  data  set. 


Such  a  rotation  matrix  can  be  defined  by  decomposing  the  covariance  matrix  T 
of  the  data  traces.  The  elements  of  the  covariance  matrix  are  given  by 


Iki  =  '^2xk(tj)xi(tj),  k,l  —  l,  (2.1) 

3= 1 


The  covariance  matrix  T  can  be  decomposed  into  its  corresponding  eigenvalues 
and  eigenvectors: 


T  =  RARr,  (2.2) 

where  A  is  a  diagonal  matrix  consisting  of  the  eigenvalues  of  T  and  R  is  the  eigenvector 
matrix  whose  columns  are  the  corresponding  eigenvectors.  Kramer  and  Mathews 
(1956)  showed  that  the  eigenvectors  define  the  principal  component  directions  and 
matrix  R7  is  exactly  a  rotation  matrix  that  decomposes  (or  rotates)  the  original 
signal  into  such  components: 

n 

fpk(t)  =  ^2  rkiXi(t)  k  =  1, ...,  n,  (2.3) 

i= 1 

and  the  original  signal  can  be  reconstructed  by: 

n 

Xi(t)  =  i  =  1, ...,  n.  (2.4) 

k=  1 

Using  a  compact  notation  with  a  data  matrix  X  whose  rows  are  the  data  traces,  the 
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decomposition  (\1/)  and  reconstruction  can  be  written  as 


\I>  =  rtx 

(2.5) 

X  =  RTq 

(2.6) 

We  note  that  the  eigenvector  matrix  R  is  a  unitary  matrix  and  satisfies: 

RRt  =  RtR  =  I. 


where  I  is  an  identity  matrix.  Thus  the  reconstruction  is  exact  and  unique,  and  the 
data  matrix  X  is  equal  to: 


X  =  X  =  RRrX.  (2.7) 

Mathematically,  the  reconstructed  signal  is  exactly  equal  to  the  original  signal  when 
all  principal  components  are  included.  Reconstructing  with  a  selected  subset  of  prin¬ 
cipal  components  ^fc(f)  would  yield  a  truncated  reconstruction  that  is  devoid  of  the 
energy  captured  in  the  discarded  components.  For  example,  incoherent  noise  often 
populates  the  last  few  principal  components  and  discarding  them  during  reconstruc¬ 
tion  would  yield  a  cleaner  signal  that  is  minimally  affected  by  the  noise.  Thus  to 
perform  simple  de-noising  of  signals,  the  reconstruction  series  is  truncated  as: 

C 

Xi(t )  =  y:rlkMt)  c  <  n,  (2.8) 

k= 1 

although  other  reconstruction  schemes  are  also  possible.  In  matrix  notation: 

X  =  rbrtx 

X  =  YX,  (2.9) 


where  B  is  a  diagonal  matrix,  with  ones  at  the  rows  corresponding  to  the  components 


used  for  reconstruction,  and  zeros  everywhere  else.  This  result  is  derived  from  the 
fact  that  PCA  measures  correlated  energy.  If  a  signal  is  random  from  trace  to  trace 
then  there  is  very  little  correlated  energy  between  the  signals,  so  the  reconstruction  of 
those  signals  falls  to  the  later  components.  For  this  reason,  truncation  of  the  rotated 
matrix  is  an  extremely  fast  process  to  remove  uncorrelated  noise  from  the  data. 
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CHAPTER  3 


ALGORITHM  DEVELOPMENT 


3.1  Introduction 

The  PCA  algorithm  for  UXO  surveys  follows  the  same  general  pattern  as  de¬ 
scribed  in  Chapter  2.  However,  significant  changes  need  to  be  made  depending  on 
whether  temporal  or  spatial  correlation  of  the  data  is  chosen.  Temporal  correlation 
requires  significant  data  regularization  before  calculation  of  the  covariance  matrix, 
while  spatial  correlation  may  require  a  distance  function  applied  to  the  correlation. 
In  addition,  with  spatial  correlation  the  covariance  matrix  can  become  extremely 
large,  precluding  a  complete  eigenvector  decomposition.  Thus  more  efficient  eigen¬ 
vector  decomposition  techniques  are  required. 

Once  decomposition  into  principal  components  is  complete,  proper  choice  of 
which  principal  components  used  in  reconstruction  is  critical.  This  choice  is  dependent 
on  the  types  of  signal  that  need  to  be  attenuated.  Since  the  orthogonal  bases  are  data- 
dependent  instead  of  user-chosen,  the  principal  components  used  in  reconstruction  are 
entirely  a  function  of  the  signal  to  be  removed,  and  are  common  throughout  all  types 
of  surveys. 

The  reconstructed  data  no  longer  has  the  same  decay  characteristics  as  before, 
so  the  same  inversion  kernels  can  not  be  used  to  invert  PCA-processed  data.  We 
therefore  developed  an  inversion  algorithm  for  both  parametric  and  generalized  inverse 
cases. 

This  chapter  describes  the  algorithm  we  developed,  from  data  normalization  to 
inversion.  In  general,  PCA  follows  the  workflow  as  described  in  figure  3.1. 
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Figure  3.1.  General  workflow  for  processing  TEM  data  with  principal  component 
analysis 


3.2  Principal  component  analysis  for  transient  electromagnetic  surveys 

In  order  to  appropriately  apply  principal  component  analysis  to  TEM  surveys, 
we  must  first  define  how  we  correlate  our  data.  We  may  either  correlate  TEM  data 
in  terms  of  space  or  in  terms  of  time.  Spatial  correlation  relates  variances  between 
decay  curves,  while  temporal  correlation  relates  variances  in  time-intersects.  For 
time  domain  surveys,  construction  of  the  covariance  matrix  must  be  modified  from 
the  standard  definition  as  described  in  equation  2.1.  The  subsequent  processing  steps 
described  earlier  may  be  directly  applied  to  this  normalized  data. 

3.2.1  Temporal  correlation 

Temporal  correlation  relates  variances  in  time-intersects.  This  method  of  orga¬ 
nization  is  extremely  useful  for  separating  signals  which  have  different  decay  charac¬ 
teristics,  such  as  UXO  and  magnetic  soils. 

While  most  multi-channel  data  can  be  decomposed  using  the  simple  covariance 
matrix  formulation  shown  in  equation  2.1,  the  particular  characteristics  of  TEM  data 
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sorted  by  time  preclude  this  construction.  Since  the  signal  decays  exponentially  with 
time,  the  data  change  in  magnitude  over  several  orders  with  a  non-linear  correspond¬ 
ing  change  in  error  level.  This  results  in  a  poorly  scaled  covariance  matrix,  which 
results  in  inaccurate  eigenvector  decomposition.  Most  of  the  covariance  is  contained 
in  the  first  few  time  gates,  overwhelming  the  later  gates.  In  addition,  for  certain  sur¬ 
veys  where  magnetic  soil  is  present  such  as  at  Kaho’olawe,  a  large  signal  (f-1  decay) 
dominates  the  data.  Figure  3.2  shows  an  example  covariance  matrix  for  an  EM-63 
survey  before  and  after  the  scaling  is  applied.  This  scaling  is  critical  for  accurate 
processing. 


In  order  to  properly  scale  the  covariance  matrix,  we  apply  normalization  to  the 
data  before  processing.  Let  X  be  defined  as  a  data  matrix  with  an  individual  row 
corresponding  to  an  individual  decay  curve  or  trace. 


Let  b  be  defined  as 

b  .  =  X,  -  (3.1) 

m 

where  m  is  the  number  of  time  gates,  Xj  is  a  row  vector  containing  an  individual 
trace,  and  n  is  the  number  of  observation  locations.  We  then  define  the  standard 
deviation  for  the  jth  time  gate  as 

/  n  \  V2 

V=(X>yJ  '  (3-2) 

We  then  normalize  the  data  by  the  first  time  gate  of  each  curve  as 


Xj 

Xi,  1 


(3.3) 


(Please  note  we  use  an  assignment  statement  rather  than  an  equality  simply  because 
we  have  run  out  of  letters.)  The  final  step  is  to  normalize  the  data,  by  the  standard 
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Figure  3.2.  Covariance  calculations  (a)  before  and  (b)  after  normalization  of  the  data. 
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deviations  calculated  in  equation  3.2: 


CTi 


(3.4) 


The  correlation  matrix  is  then  calculated  from  this  normalized  data.  As  described 
in  detail  in  Chapter  2,  the  data  undergo  an  orthogonal  transformation  to  a  basis 
spanned  by  the  eigenvectors  of  this  new  covariance  matrix.  Principal  components 
are  chosen  to  be  removed  (the  choice  of  which  are  described  in  sections  3.3  and  3.4), 
and  the  data  are  rotated  back,  leaving  a  dataset  with  unwanted  features  removed  or 
attenuated.  It  is  important  to  note  that  at  the  end  of  PCA  processing,  the  data  must 
be  un-normalized  by  applying  the  above  procedure  in  reverse  order  to  the  processed 
data. 


3.3  Removal  of  uncorrelated  noise 

The  removal  of  uncorrelated  noise  generally  involves  the  removal  of  the  last  few 
principal  comp onents-t hose  corresponding  to  the  smallest  eigenvalues.  The  choice 
of  how  many  principal  components  to  remove  can  be  made  by  following  a  general 
guideline  based  on  the  estimated  noise  threshold.  The  principal  components  are 
removed  until  the  sum  of  their  associated  eigenvalues  relative  to  the  total  is  equal  to 
the  estimated  noise  level.  For  example,  if  we  estimate  5%  random  noise  in  our  survey 
data,  then  we  remove  principal  components  starting  with  the  last  one  until  the  sum 
of  the  associated  eigenvalues  divided  by  the  total  sum  of  the  eigenvalues  equals  0.05. 

These  principal  components  are  removed  by  applying  a  truncation  matrix  to  the 
rotated  data,  as  described  in  equation  2.9.  The  truncation  matrix  starts  as  an  identity 
matrix,  where  elements  on  the  diagonal  are  replaced  with  zeros  corresponding  to  the 
principal  components  to  be  removed. 
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3.4  Removal  of  correlated  noise 


The  process  of  removal  of  correlated  noise  is  dependent  on  the  type  of  noise  to 
be  removed.  This  type  of  noise  can  be  classified  into  two  types,  each  having  a  distinct 
and  well-defined  method  for  removal:  static  and  variable. 

Static  noise  includes  sources  such  as  magnetic  soils,  which  exhibit  decay  pro¬ 
portional  to  t~1,  and  sensor  orientation,  which  behaves  like  a  static  shift  for  small 
perturbations  of  the  orientation  of  receiver /transmitter  coils.  These  sources  of  noise 
tend  to  effect  large  areas  of  the  survey  grid  and  have  a  temporally-invariant  signa¬ 
ture.  These  noise  sources  map  into  the  first  principal  component,  thus  removal  of  the 
first  principal  component  will  attenuate  these  signals.  Removal  of  the  first  principal 
component  will  also  remove  a  portion  of  the  signal  due  to  UXO,  but  will  leave  the 
diagnostic  decay  characteristics  intact  (Figures  4.8  and  4.10).  Section  3.8  describes 
how  to  interpret  the  data  after  processing  while  incorporating  this  effect. 

Variable  noise  includes  signals  from  telluric  currents,  RF  interference,  and  similar 
signals.  These  are  more  temporally  and  spatially  variable,  and  will  thus  map  into 
principal  components  smaller  than  those  containing  the  UXO.  Therefore,  in  surveys 
where  noise  sources  such  as  these  are  known  to  be  present,  only  the  first  two  principal 
components  can  generally  reconstruct  the  signal  due  to  UXO,  while  the  signal  due  to 
this  noise  is  attenuated. 

3.5  Spatial  correlation 

Spatial  covariances  relate  changes  in  anomaly  geometry  across  a  survey  rather 
than  changes  in  decay  characteristics.  While  this  method  of  organization  is  less 
useful  for  UXO  applications  than  temporal  correlation,  the  method  can  still  be  useful 
in  applications  for  removal  of  large  geologic  features. 

Initial  construction  of  the  covariance  matrix  requires  no  special  normalization 
of  the  data  as  the  variation  in  measured  magnitude  across  a  single  time-intersect  is 
much  smaller  than  across  a  decay  curve.  However,  choice  of  principal  components  to 
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remove  is  no  longer  clear  and  requires  significant  user  input.  How  a  particular  UXO 
item  maps  into  the  principal  component  space  is  highly  dependent  on  the  spatial 
characteristics  of  the  geology-much  more  so  than  in  a  temporal  correlation.  Thus 
trial  and  experience  are  critical  for  proper  choice  of  principal  components  for  recon¬ 
struction.  As  a  consequence  of  this  dependence,  UXO  can  map  into  a  much  wider 
variety  of  principal  components.  Therefore  calculating  the  first  few  eigenvectors  to 
save  computation  time  is  not  an  appropriate  approach  to  the  problem.  Since  the  co- 
variance  matrix  in  this  construction  has  n2  elements,  where  n  is  the  number  of  data 
traces,  the  matrices  can  be  extremely  large. 

For  these  reasons,  we  do  not  recommend  the  use  of  spatially-correlated  covari¬ 
ances  as  a  component  of  a  standard  UXO  workflow.  However,  cases  may  exist  where 
strong  spatially- variant  features  have  the  same  decay  characteristics  as  UXO,  and 
the  spatial  correlation  must  be  taken  advantage  of.  Section  3.6  discusses  numerical 
techniques  to  address  these  issues. 


3.6  Large-scale  decompositions 


When  correlating  data  spatially,  the  covariance  matrices  can  become  extremely 
large.  Thus  efficient  methods  of  eigenvector  decomposition  must  be  used.  Two  basic 
approaches  can  be  taken:  decimating  the  covariance  matrix  to  make  it  sparse,  or 
utilizing  efficient  eigenvector  decomposition  algorithms. 

Decimating  the  covariance  matrix  involves  assuming  that  certain  data  locations 
must  have  zero  correlation,  usually  due  to  distance.  Applying  a  distance  weighting 
function  to  the  covariance  calculation,  whether  smooth  or  binary,  can  result  in  a 
sparse  matrix.  As  an  example  of  a  smooth  weighting  function,  we  apply  a  distance 
weight  to  the  covariance  calculation: 


'y ki  ^  ^  %kj%ij  ( 1 


3= 1 


ay/(yk  ~  Vi )2  +  (zk  ~  zi)2 
s/max(yi  -  y0 )2  +  rna x(zt  -  z0 )2 


i,k,l  —  1, ...,  n  (3.5) 
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where  y ,  z  are  Cartesian  coordinates  of  each  data  location,  and  a  is  some  constant 
that  controls  the  radius  of  inclusion.  By  choosing  an  appropriate  algorithm  that 
utilizes  this  sparsity,  the  computation  time  and  memory  requirements  can  be  greatly 
reduced. 

Because  of  the  size  and  potentially  sparse  nature  of  the  spatial  covariance  matri¬ 
ces,  iterative  methods  for  eigenvector  decompositions  are  most  appropriate  for  PCA. 
There  are  many  to  choose  from,  each  with  advantages  and  disadvantages,  and  no  one 
decomposition  is  perfect  in  all  cases.  However,  we  briefly  review  the  most  applicable 
iterative  methods  to  UXO  here.  For  further  details  on  each  of  these  algorithms,  we 
refer  the  reader  to  Blackford  (2000). 


3.6.1  Iterative  eigenvalue  decompositions 

Lanczos  method 

The  Lanczos  method  is  one  of  the  most  widely  used  algorithms  for  large-scale 
eigenvector  decomposition.  By  building  up  Ritz  approximations  to  the  eigenvalues 
of  the  matrix,  the  Lanczos  method  produces  several  eigenvalues  from  once  sequence 
of  vectors  that  converge  quickly.  Unfortunately,  even  mildly  ill-conditioned  matrices 
produce  a  system  which  destroys  orthogonality  quickly,  requiring  a  reorthogonaliza- 
tion  process  (i.e.  Gram-Schmidt). 


Power  method 

The  power  method  is  the  most  simple  of  the  methods  to  implement.  However, 
the  process,  like  many  of  the  other  methods,  begins  at  one  extreme  eigenvalue  and 
calculates  each  subsequent  eigenvalue-finding  a  subset  of  interior  eigenvalues  requires 
a  shift-and-invert  modification. 
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Subspace  iteration  method 


This  method  is  extremely  common  in  engineering  applications,  and  may  be  the 
most  readily  applied  process  to  PCA  due  to  its  prevalence  in  the  community.  Much 
like  the  power  method,  however,  the  subspace  iteration  method  directly  produces 
only  leading  eigenvalues  and  eigenvectors. 

Jacobi-Davidson  method 

The  Jacobi-Davidson  method  can  directly  produce  a  cluster  of  interior  eigenval¬ 
ues  without  calculating  the  entire  set  by  using  preconditioners,  making  it  extremely 
useful  in  spatially-correlated  matrices. 

3.7  Anomaly  selection 

Once  PCA  has  been  performed  on  a  dataset,  the  data  are  ready  for  processing 
with  automated  anomaly  picking  routines  such  as  with  UXOLab,  or  manual  interpre¬ 
tation.  Each  of  the  chosen  anomalies  can  then  be  further  processed  with  numerical 
inversion  techniques  or  other  discrimination  algorithms. 

3.8  Inversion  of  processed  data 

While  inversion  of  electromagnetic  data  is  a  nonlinear  process,  insight  into  the 
requirements  for  inversion  of  PCA-processed  data  can  be  gleaned  from  studying  a 
linear  as  well  as  nonlinear  case.  Linear  inversion  of  geophysical  data  seeks  to  invert 
the  forward  mapping  operator  (or  sensitivity  matrix),  G  that  operates  on  the  model 
to  produce  measured  data.  This  inverted  matrix  may  then  be  applied  to  the  data 
to  map  them  back  into  the  model  space.  In  general,  noise  must  be  accounted  for 
and  the  problem  is  ill-posed.  Thus  further  information,  such  as  constraints  on  the 
size  or  structure  of  the  model,  reference  models,  data  weighting,  and  many  others  is 
required  for  a  stable  inversion  (Oldenburg  and  Li,  2005).  Allowing  a  trade  off  between 
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the  model  norm  ((f>m)  and  the  data  misfit  (4>d),  leads  to  the  global  objective  function 
which  is  minimised  to  arrive  at  the  Tikhonov  solution: 


0  —  +  /34>m- 


(3.6) 


Inverting  for  the  smallest  model  with  an  L2  data  misfit  measure,  the  minimization  of 
this  solution  expands  to: 

(GTG  +  /?I)m  =  GTd  (3.7) 

and  m  can  be  solved  for  either  directly  or  iteratively  using  the  conjugate  gradient  or 
similar  solution.  The  Tikhonov  or  regularization  parameter,  (3,  is  either  chosen  such 
that  an  optimal  data  misfit  is  reached,  by  using  the  L-curve  criterion  (Hansen,  1992), 
or  a  generalized  cross  validation  (GCV)  approach. 

The  forward  mapping  operator  described  here  encompasses  the  physics  and  ge¬ 
ometry  of  the  geophysical  problem.  Once  data  has  been  processed  with  PCA  for 
signal  isolation  or  de-noising,  this  forward  mapping  operator  no  longer  accurately 
maps  the  model  to  this  new  rotated  data  as  the  data  have  been  rotated  to  a  new 
basis  in  Rn.  In  order  to  accomplish  this  mapping  (and  thus  inversion),  the  forward 
mapping  operator  must  also  be  rotated  such  that  it  maps  from  the  model  space  to 
this  new  data  space. 


We  apply  the  rotation  matrix  to  each  column  of  the  sensitivity  matrix,  G,  in¬ 
dividually  to  match  the  processed  data,  dr0i.  This  is  equivalent  to  multiplying  a 
new  matrix  Yg  to  the  sensitivity  matrix,  where  Yg  is  a  block-sparse  matrix  with  the 
original  rotation  matrix  Y  on  the  diagonal  as: 


Y 


a 


Y  0  ...  0 

0  Y 


0 


Y 


(3.8) 
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Thus  Equation  3.7  becomes 


((Y9G)t(Y9G)  +  /3I)m  =  (Y9G)Tdro,  (3.9) 

and  can  be  solved  with  the  same  numerical  methods  as  before.  It  is  important  to 
note  that  Y  and  Yg  are  both  positive-semidefinite  since  they  are  calculated  from  the 
eigenvectors  of  the  covariance  matrix. 

The  general  solution  to  a  non-linear  inverse  problem  is  solved  by  developing  a 
locally  linear  system  and  minimising  the  objective  function  described  in  Equation  3.6. 
We  then  update  our  sensitivity  matrix  and  continue  in  an  iterative  process  until  the 
global  objective  function  is  minimised  via  a  Gauss-Newton  or  similar  method.  In  the 
case  where  no  rotation  on  the  data  has  been  performed,  we  may  develop  our  global 
objective  function  with  a  data  misfit  and  model  norm  to  be  minimised  which  expands 
to: 

0(m)  =  ||Wd(do6s  -  E[m])||2  +  /3||Wm(m  -  mre/)||2  (3.10) 

where  m  is  the  vector  of  parameters,  dcbs  is  the  recorded  data  with  a  corresponding 
weighting  matrix  W^,  F  is  the  forward  mapping  operator  (that  depends  on  m),  (3 
is  a  Tikhonov  or  regularization  parameter,  m ref  is  a  reference  model,  and  Wm  a 
corresponding  weighting  matrix. 

In  general,  the  local  linear  system  can  be  obtained  through  a  perturbation  ap¬ 
proach  using  the  Gauss-Newton  method  with  the  following  form: 

(JTWJW,J  +  pW^Wm)  5m  = 

JTWjW,  ( dobs  -  F  [m^])  -  pWlWm  (mW  -  mre/)  (3.11) 

where  J  is  a  Jacobian  matrix.  In  order  to  generalize  the  forward  mapping  such  that 
the  rotation  and  possible  truncation  of  the  data  is  included,  we  apply  the  same  rota¬ 
tion  operator  to  the  data  predicted  by  the  forward  operator  F,  just  as  in  the  linear 
case.  This  leads  to  an  altered  equation  3.11: 
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(3.12) 


(jLwJWdJroi  +  /3W^Wm)  5m  = 

JrotWjWd  (drot  -  Y,F  [mW])  -  /?W£Wm  (mW  -  mre/) 

where  drot  is  the  rotated  data  and  Yg  is  the  same  matrix  defined  in  Equation  3.8. 
The  inversion  can  be  computed  in  the  same  way  as  before.  Since  the  Jacobian  matrix 
is  calculated  from  the  predicted  data,  the  rotation  matrix  will  be  combined  into  the 
Jacobian,  Jrot,  with  each  subsequent  step,  so  explicit  rotation  of  the  Jacobian  is 
unnecessary. 

Alteration  of  the  original  data  by  a  transform  also  modifies  the  statistical  dis¬ 
tribution  of  the  noise.  After  PCA,  ideally  there  is  no  random  noise  and  all  signal 
is  a  result  of  coherent  sources  (whether  geologic  or  otherwise).  This  would  imply 
that  in  an  ideal  separation  case,  the  Tikhonov  parameter  f3  would  approach  zero  to 
maximally  weight  a  data  misfit  of  zero.  Clearly  this  case  is  not  achieved  for  several 
reasons:  data  errors  are  not  always  zero  bias  with  a  Gaussian  distribution  (or  any 
other  uncorrelated  distribution),  PCA  has  intrinsic  numerical  error  in  the  calculation 
of  eigenvectors,  there  is  strongly  correlated  noise  associated  with  most  surveys,  and 
original  statistical  estimates  for  error  removal  may  have  been  incorrect  in  the  Erst 
place  (i.e.  which  principal  components  to  use).  Therefore  the  data  misfit  is  guaran¬ 
teed  to  not  have  a  \2  distribution.  So  while  the  value  of  the  data  misfit  is  defined, 
the  appropriate  value  of  data  misfit  for  an  optimal  solution  is  not.  However,  this 
altered  data  misfit  can  be  treated  while  choosing  f3  with  the  L-curve  criterion  or  with 
a  generalized  cross  validation  approach-both  work  regardless  of  the  rotation  applied. 

The  same  approach  to  inverting  rotated  data  in  a  generalized  inversion  applies 
equally  to  a  parametric  inversion.  In  the  parametric  case,  the  global  objective  function 
does  not  contain  a  model  objective  (as  the  model  is  defined  explicitly  already),  and 
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1 

Perform  any  initial  filtering  (such  as  median  filter¬ 
ing),  if  necessary 

2 

Normalize  data  (3.2.1) 

3 

Data  decomposition  (2.1) 

4 

Principal  component  truncation  (2.1,  3.3,  and  3.4) 

5 

Reconstruction  (2.1) 

6 

Remove  data  normalization  (3.2.1) 

7 

Application  of  anomaly  automatic  picking  rou¬ 
tines,  such  as  in  UXOLab  (3.7) 

8 

Numerical  inversion  of  individual  anomalies  (3.8) 

Table  3.1.  Recommended  PCA  algorithm.  Numbers  in  parentheses  correspond  to 
sections  with  detailed  descriptions. 


so  a  Gauss-Newton  solution  to  a  parametric  case  with  rotated  data  becomes 

(JTWJW,J)  5m  =  JTWjWd  ( dobs  -  Y gF  [mW])  . 

As  in  the  full  generalized  inversion,  the  rotation  matrix  should  not  be  explicitly 
applied  to  the  Jacobian. 

3.9  Algorithm  summary 

In  this  chapter  we  have  introduced  a  PCA  algorithm  for  processing  TEM  data  to 
attenuate  both  uncorrelated  and  correlated  noise,  as  well  as  developed  an  inversion 
algorithm  for  interpretation  of  individual  UXO  anomalies.  Table  3.1  shows  a  skeleton 
algorithm  for  the  processing  and  interpretation  of  TEM  data  with  PCA.  In  the  next 
chapter,  we  show  examples  of  the  successful  application  of  PCA  to  TEM  data. 
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CHAPTER  4 


RESULTS 


In  this  chapter  we  present  three  examples  of  the  successful  application  of  PCA  to 
unexploded  ordnance  detection  and  discrimination.  These  examples  show  the  varied 
applications  and  contributions  of  PCA  to  detection  and  discrimination  of  UXO,  as 
well  as  the  ease  of  insertion  into  a  standard  workflow. 

4.1  Uncorrelated  noise  removal 

To  show  the  de-noising  capability  of  PCA,  we  present  an  example  dataset  from 
Kaho’olawe,  Hawaii,  USA.  The  TEM  data  (collected  with  a  Geonics  EM63)  come 
from  a  30  metre  by  60  metre  grid  set  on  extremely  magnetic  soils  (magnetite,  titano- 
magnetite,  and  ilmcnite)  due  to  volcanism.  Weathering  of  the  volcanics  in  the  tropical 
climate  led  to  soils  exhibiting  a  high  degree  of  frequency  dependence  in  magnetic  sus¬ 
ceptibility.  The  spatially  variable,  frequency-dependent  susceptibility  produces  strong 
1/t  decay  that  masks  the  UXO  response  in  EM63  data.  The  soil  response  is  over  400 
mV  in  the  first  gate  across  the  survey  area  (Figure  4.1).  In  addition,  the  instrument 
movement,  ambient  interferences,  and  other  sources  lead  to  the  ubiquitous  noise  in 
late  time  gates. 

To  assist  with  the  separation  of  the  magnetic  soil  response  from  unexploded 
ordnance,  we  must  first  attenuate  the  noise  in  the  late  time  gates.  The  strong  soil 
response  masks  much  of  the  signal  due  to  UXO  before  the  last  few  time  gates,  so 
improvement  in  signal-to-noise  is  paramount  in  interpretation. 

The  data  (consisting  of  23000  decay  curves  of  25  gates  each)  were  decomposed 
into  constituent  principal  components  and  re-composed  using  a  subset  of  those  com¬ 
ponents.  Figure  4.2(a)  shows  a  set  of  the  original  data  that  include  only  geologic 
response  as  well  as  data  that  have  signal  from  UXO  imprinted  on  the  decay  curves. 
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Figure  4.1.  First  time  gate  of  EM63  data  collected  over  the  Kaho’olawe  Navy  QA 
grid.  Note  the  large  geologic  response  due  to  viscous  remnant  magnetization  (VRM). 
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Figure  4.2.  (a)  65  example  decay  curves  from  Kaho’olawe.  Some  traces  include 

responses  from  unexploded  ordnance;  these  traces  deviate  from  the  expected  t 
decay,  (b)  Geologic  response  constructed  with  only  the  first  principal  component.  The 
curve  exhibits  the  expected  inverse  power  law  decay  due  to  horizontal  layering  and 
magnetic  soils,  (c)  TEM  decay  curves  constructed  with  the  first  and  second  principal 
components.  The  (approximate)  superposition  of  UXO  signal  on  the  geologic  signal 
is  clear,  with  significant  reduction  of  random  noise  in  later  channels.  (Negative  values 
not  displayed) 
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Magnetic  soil  at  Kaho’olawe  (and  in  general)  exhibits  a  t~l  decay  in  TEM  surveys 
(Pasion  et  ah,  2002).  This  signal  is  clearly  present  in  the  raw  data  (Figure  4.2(a)). 
Reconstruction  with  only  the  first  principal  component  reveals  this  decay  due  to  the 
magnetic  soil  (Figure  4.2(b)).  Traces  with  UXO  and  those  without  UXO  are  indistin¬ 
guishable.  However,  when  the  second  principal  component  is  added  (Figure  4.2(c)), 
the  signal  from  the  UXO  (deviation  from  t~l  decay  in  later  time  gates)  is  clear. 
The  random,  uncorrelated  noise  that  contaminated  these  later  time  gates  is  severely 
reduced  without  attacking  the  amplitudes  of  the  decay  curves,  preparing  them  for 
numerical  discrimination  analysis  with  other  methods. 

4.2  Inversion  of  processed  data 
4.2.1  Linear  inversion 

As  an  example  of  the  need  to  include  the  rotation  matrix  in  inversion,  we  present 
the  following  linear  example.  The  model,  m,  consisting  of  one-thousand  elements,  is 
represented  by  a  sine  function  of  the  form 

m(x)  =  1  —  .5[cos(27nr)  +  sin(27nr)] 

with  10  channels  of  data  calculated  at  each  of  11  data  locations  simulated  (110  total 
data).  The  channels  all  utilize  exponentials  as  kernels  with  different  decay  parame¬ 
ters  and  are  functions  of  current  and  adjacent  observation  locations  with  5%  Gaussian 
noise.  Thus  we  have  a  linear,  underdetermined  system  to  invert  for  the  model  pa¬ 
rameters. 

For  the  first  case  where  no  principal  component  analysis  was  applied,  we  in¬ 
verted  the  data  for  the  smallest  model.  The  Tikhonov  parameter,  /3,  was  chosen 
through  the  L-curve  criterion.  Because  our  chosen  kernels  are  sensitive  to  noise,  the 
smallest-model  solution  yields  a  poor  recovered  model  (Figure  4.3(a)).  Even  though 
the  optimal  data  misfit  is  reached,  the  structure  of  the  recovered  model  is  poorly  con- 


strained  at  depth.  The  imprinted  noise  on  the  underdetermined  problem  is  preventing 
a  good  result,  despite  an  excellent  L-curve  (Figure  4.4(a)). 

As  a  de-noising  tool,  principal  component  analysis  is  a  logical  choice  to  reduce 
the  uncorrelated  noise  present  in  the  data.  Therefore,  we  applied  a  blind  PCA  to  the 
data,  where  the  amount  of  energy  contained  in  the  noise  was  estimated  and  latter 
principal  components  were  removed  to  reach  5%  noise,  nulling  all  but  the  first  two 
principal  components.  The  data  were  then  rotated  back  and  inverted  again  with  the 
same  sensitivity  matrix  as  before  to  yield  the  model  in  Figure  4.3(b).  Though  the 
shallow  structure  has  improved,  it  still  does  not  accurately  represent  the  true  model. 
In  addition,  choice  of  regularization  parameter  was  more  ambiguous  than  the  previous 
case,  as  the  L-curve  contained  two  areas  of  high  curvature  (Figure  4.4(b)). 

To  improve  the  recovered  model,  we  applied  the  rotation  matrix  to  the  sensitivity 
matrix  as  described  in  this  paper.  With  the  sensitivity  matrix  consistent  with  the 
data  space,  the  inversion  produced  a  much  better  result  (Figure  4.3(c)).  Moreover, 
the  optimal  f3  term  was  easy  to  choose  from  the  L-curve  (Figure  4.4(c)). 

4.2.2  Non-linear  inversion 

We  present  synthetic  data  from  an  unexploded  ordnance  survey.  We  simulate 
data  from  a  25  channel  transient  electromagnetic  survey  (such  as  with  a  Geonics  EM- 
63,  for  example)  over  a  half-space  with  a  75mm  mortar  round  buried  at  1  m  depth. 
We  use  the  Pasion-Oldenburg  parametric  model  that  uses  two  dipolar  responses  to 
simulate  the  signal  from  the  mortar  round.  Each  dipole  decays  as: 

k{t  +  a)~^e~th  (4.1) 

where  the  parameters  k ,  a,  f3,  and  7  depend  on  the  conductivity,  permeability,  size, 
and  shape  of  the  object.  The  ratios  of  certain  values  can  then  be  used  as  diagnostic 
indicators  for  unexploded  ordnance  (LIXO)  candidates. 

We  parametrically  invert  the  synthetic  dataset  for  the  LIXO  parameters  with 
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Figure  4.3.  Original  and  recovered  model  (a)  through  smallest  model  solution,  (b) 
after  PCA  on  data  only,  and  (c)  after  PCA  on  data  and  kernels. 
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Figure  4.4.  L-curve  for  (a)  original  inversion,  (b)  inversion  with  rotated  data,  and  (c) 
inversion  with  rotated  data  and  rotated  sensitivity  matrix. 
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Table  4.1.  Nonlinear  inversion  results.  Only  by  applying  PCA  to  both  the  data  and 
the  kernels  were  we  able  to  properly  recover  diagnostic  parameters. 


the  intention  of  recovering  the  proper  diagnostic  ratios.  The  objective  function  is 
constructed  in  the  same  manner  as  before,  save  that  there  is  no  model  weighting 
(since  this  is  a  parametric  inversion).  Thus  the  linearized  system  to  solve  becomes: 

( JTWjW(/J)hm  =  JTWjW,(dobs  -  F  [mW] )  .  (4.2) 

We  solve  this  system  for  noisy  data,  PCA-processed  data,  and  PCA-processed 
data  with  the  rotation  matrix  incorporated  into  the  inversion.  Table  1  shows  the 
results.  Inverting  the  data  with  no  processing  fails  to  recover  the  k  ratio,  which 
contains  shape  information.  Processing  the  data  with  PCA  but  not  incorporating  the 
rotation  matrix  into  the  inversion  results  in  incorrect  recovery  of  all  shape  information. 
Only  by  incorporating  the  rotation  matrix  do  we  successfully  identify  a  potential  UXO 
target. 

4.3  Correlated  noise  removal 

Variations  in  the  orientation  of  the  TEM  transmitter  and  receiver  cart  due  to 
microtopography  can  produce  significant  signal  in  the  data.  Figure  4.5  shows  the  first 
time  gate  recorded  in  a  synthetic  EM-63  survey.  This  survey  used  real  orientation 
vectors  from  a  UXO  survey  and  simulated  the  TEM  response  of  the  orientation  errors 
and  40  mm  mortar  shells.  Large,  linear  trends  can  be  seen  in  the  data  where  an 
approximate  static  shift  in  amplitude  across  all  time  gates  has  been  introduced. 

This  type  of  correlated  noise  often  has  a  large  effect  on  autopicking  routines. 
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Figure  4.5.  Effect  of  cart  orientation  error  in  a  UXO  survey.  Orientation  error 
produces  an  approximate  50  mV  response  in  this  example,  visible  in  the  approximate 
east/west  lineations. 


The  characteristics  of  correlated  noise  sources  are  often  similar  enough  to  those  of 
UXO  that  autopicking  routines  cannot  differentiate  them.  Moreover,  the  frequency 
contents  of  the  signals  overlap  such  that  frequency  filtering  is  difficult  or  impossible. 
Figure  4.6  shows  the  radially-averaged  power  spectra  of  the  constituent  components  in 
the  dataset  shown  in  Figure  4.5.  The  power  spectra  of  signal  due  to  cart  orientation 
error  concentrates  power  in  the  same  bands  as  the  signal  due  to  UXO.  Therefore 
separation  using  simple  frequency  filtering  is  impossible  (Figure  4.7).  Processing  the 
dataset  in  Figure  4.5  with  UXOLab’s  (UXOLab,  2009)  autopicking  routines  leads  to 
60  targets  (after  clustering  the  picks).  Figure  4.8  shows  the  picked  anomalies  using 
the  routines  provided  in  UXOLab.  The  effect  of  the  cart  orientation  on  the  picking 
routine  is  clear-large  numbers  of  false  positives  result. 

Effects  due  to  cart  motion  effectively  result  in  an  approximate  static  shift  in  the 
data.  Therefore  the  orientation  error  effects  each  time  gate  equally.  Thus  to  remove 
this  effect,  the  first  principal  component  was  removed.  Figure  4.10  shows  the  dataset 
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Figure  4.6.  Radially  averaged  power  spectrum  of  the  original  data  set,  orientation 
error,  and  UXO  signal.  The  difference  between  the  power  spectra  of  the  orientation 
error  and  the  UXO  signal  indicates  that  they  are  inseparable  in  the  frequency  domain. 


with  the  cart  orientation  error  removed.  Though  the  amplitude  of  the  UXO  anomalies 
have  been  reduced,  the  shape  has  been  preserved. 

This  dataset  was  then  processed  in  UXOLab  for  automatic  UXO  anomaly  pick¬ 
ing.  Figure  4.9  shows  two  example  lines  from  one  time  gate  in  this  dataset.  Fig¬ 
ure  4.9(a)  shows  a  line  that  passes  through  signals  due  to  both  UXO  and  cart  orien¬ 
tation  error.  For  autopicking  purposes,  a  threshold  of  50  mV  was  chosen  based  on  the 
signals  present  in  this  line.  This  resulted  in  the  picks  previously  shown  in  Figure  4.8. 
After  principal  component  analysis,  a  threshold  of  8  mV  was  chosen  (Figure  4.9(b)). 
Though  visually  this  line  does  not  appear  to  be  improved,  the  new  threshold  dramati¬ 
cally  improved  the  autopicking  ability.  Figures  4.9(c)  and  4.9(d)  show  a  line  that  does 
not  pass  through  a  UXO  target,  but  does  pass  through  signal  due  to  cart  orientation 
error.  The  signal  due  to  orientation  error  is  completely  eliminated. 

While  the  original  dataset  had  62  picks  after  clustering  (Figure  4.8),  the  pro- 
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Figure  4.7.  Three  different  examples  of  Butterworth  filtering.  The  filtering  is  unable 
to  separate  the  cart  orientation  error  and  UXO  for  a  variety  of  cutoff  wavenumbers. 
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Figure  4.8.  UXOLab  target  picks  in  a  dataset  contaminated  with  orientation  error. 
Note  the  large  number  of  picks  corresponding  to  the  cart  orientation  error  (linear 
bands  of  picks).  The  dashed  line  corresponds  to  transects  shown  in  Figure  4.9. 
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Figure  4.9.  Example  lines  from  one  time  gate  before  and  after  PCA.  Units  in  mV. 
(a)  Line  including  both  UXO  anomalies  as  well  as  cart  orientation  error,  (b)  Same 
line  after  PCA.  Choosing  a  threshold  from  this  line  (8  mV)  produced  better  anomaly 
picking  results  than  from  the  non-PCA  line.  The  numbers  assigned  to  each  peak 
are  assigned  by  the  autopicking  routine,  and  thus  change  before  and  after  PCA. 
(c)  Example  line  that  contains  no  UXO  anomalies.  Note  the  large  anomaly  due  to 
orientation  error,  (d)  Same  line  after  PCA.  No  appreciable  structure  is  present,  and 
the  magnitude  is  well  below  that  of  the  UXO  anomalies. 
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Figure  4.10.  UXOLab  anomaly  picks  after  processing  with  PCA.  After  clustering, 
there  was  one  false  anomaly  and  no  missed  targets. 


cessed  dataset  produced  32  picks,  only  one  of  which  was  false.  There  were  no  missed 
targets  (Figure  4.10).  This  dataset  is  now  ready  for  further  processing  via  inversion 
of  individual  chosen  anomalies  or  other  methods.  With  standard  techniques,  these 
locations  can  be  used  on  the  original  data  for  further  processing,  ffowever,  in  order 
to  incorporate  PCA,  UXOLab  and  other  UXO  toolboxes  should  be  further  developed 
to  add  rotation  into  their  subsequent  processing  and  inversion. 
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CHAPTER  5 


SUMMARY 


This  is  the  final  report  for  the  SERDP  project  MM-1640  and  it  covers  the  re¬ 
search  results  accomplished  since  the  inception  of  the  project  in  2007.  The  basic 
premise  of  this  project  is  the  theoretical  understanding  and  algorithm  development 
of  principal  component  analysis  as  a  denoising  and  signal-separation  tool  for  transient 
electromagnetic  (TEM)  data  in  unexploded  ordnance  (UXO)  applications.  There  is 
an  express  need  for  techniques  to  reduce  the  presence  of  random  noise  in  TEM  data 
as  well  as  reduce  the  influence  of  correlated  noise  due  to  a  wide  variety  of  sources 
on  automatic  anomaly-picking  routines  for  more  accurate  detection  with  fewer  false 
anomalies.  We  have  developed  an  algorithm  and  workflow  for  the  processing  and 
inversion  of  TEM  data  that  attenuates  signal  from  undesired  sources  and  accurately 
inverts  TEM  data  for  diagnostic  UXO  parameters. 

We  separated  the  research  in  this  project  into  aspects  of  PCA  that  address 
uncorrelated  signals  and  aspects  that  address  correlated  but  undesired  signals,  as 
well  as  interpretation  of  processed  data.  Specifically,  these  objectives  included: 

•  Develop  a  PCA  algorithm  tailored  to  TEM  data  from  UXO  surveys,  and 

•  Produce  a  stable  and  automated  algorithm  for  removal  of  uncorrelated  noise. 

•  Study  the  physical  connection  between  signals  due  to  different  sources  and  the 
components  recovered  from  PCA 

•  Determine  the  feasibility  of  using  PCA  to  reduce  TEM  data  in  UXO  applications 
to  the  signal  exclusively  produced  by  UXO  and  UXO-like  anomalies,  and 

•  Apply  the  developed  algorithm  to  data  examples  and  assess  the  effectiveness 
through  discrimination/classification  tests. 
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•  Develop  a  method  for  parametric  inversion  of  processed  data  for  recovery  of 
diagnostic  parameters  in  UXO  applications. 

We  have  successfully  accomplished  these  research  goals  over  the  course  of  the 
project. 

Develop  a  PCA  algorithm  tailored  to  TEM  data  from  UXO  surveys 

We  developed  an  algorithm  based  on  the  Karhunen-Loeve  transform  tailored  to 
UXO  surveys  by  producing  a  temporally-correlated  covariance  matrix  with  proper 
data  scaling.  We  have  shown  that  this  data  organization  and  data  scaling  produces 
an  effective  algorithm  for  processing  of  UXO  data.  Spatially-correlated  covariances 
investigate  different  characteristics  of  TEM  data  and  are  not  immediately  well-suited 
for  UXO  processing.  In  addition,  spatial  correlations  require  intensive  numerical 
processing  and  significant  user  intervention  in  the  process. 

Produce  a  stable  and  automated  algorithm  for  removal  of  uncorrelated 
noise 

The  PCA  algorithm  we  developed  reliably  and  automatically  removes  uncorre¬ 
lated  noise.  Using  an  estimate  of  the  error  level  in  the  data,  the  PCA  algorithm 
automatically  removes  the  proper  principal  components  that  contain  the  signal  clue 
to  this  uncorrelated  noise.  This  is  accomplished  with  a  truncation  matrix  applied  to 
the  inverse  rotation  matrix  that  is  defined  by  the  number  of  eigenvalues  required  to 
reach  the  error  level  estimate.  The  processed  data  contain  unaltered  signal  clue  to 
UXO  with  the  random  noise  strongly  attenuated. 

Study  the  physical  connection  between  signals  due  to  different  sources  and 
the  components  recovered  from  PCA 

Undesired  signals  in  TEM  surveys  map  into  different  principal  components  as 
a  function  of  their  decay  characteristics.  Both  numerical  simulation  and  field  data 
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confirm  that  sources  which  contain  static  shifts,  such  as  sensor  orientation  error  and 
magnetic  soils,  map  into  the  first  principal  component.  These  static  sources  tend  to 
overlap  with  some  of  the  signal  produced  by  UXO;  however,  even  with  truncation  of 
the  first  principal  component,  diagnostic  decay  and  significant  anomalies  remain  due 
to  the  UXO,  and  do  not  degrade  performance  of  automated  anomaly  picking  routines, 
such  as  those  available  in  UXOLab.  Sources  which  produce  temporally- variant  signals 
have  much  more  variety  in  how  they  map  into  principal  components,  but  tend  to  be 
well-separated  from  UXO  and  map  into  the  3rd  and  higher  principal  component. 


Determine  the  feasibility  of  using  PCA  to  reduce  TEM  data  in  UXO  appli¬ 
cations  to  the  signal  exclusively  produced  by  UXO  and  UXO-like  anomalies 

We  have  determined  through  testing  and  experimentation  of  the  algorithm  on 
synthetic  and  held  data  that  PCA  can  separate  signals  due  to  correlated  noise  from 
signals  produced  by  UXO  and  UXO-like  anomalies.  We  have  also  determined  that  the 
temporally-correlated  PCA  algorithm  developed  here  is  not  appropriate  for  clutter 
analysis  of  individual  anomalies  when  the  entire  dataset  is  processed  at  once.  How¬ 
ever,  current  research  in  SERDP  and  other  projects  suggests  that  various  types  of 
component  analysis  are  effective  in  small,  windowed  areas  for  separating  signals  due 
to  closely  spaced  targets. 


Apply  the  developed  algorithm  to  data  examples  and  assess  the  effective¬ 
ness  through  discrimination/classification  tests 

We  have  successfully  applied  PCA  to  several  synthetic  and  held  datasets.  In 
all  cases,  the  signal-to-noise  ratio  was  increased,  the  ability  to  detect  anomalies  in 
the  presence  of  correlated  noise  was  improved,  and  numerical  inversion  of  individual 
anomalies  of  processed  data  was  made  possible. 


41 


Develop  a  method  for  parametric  inversion  of  processed  data  for  recovery 
of  diagnostic  parameters  in  UXO  applications 

Using  the  Pasion-Oldenburg  model  of  the  TEM  response  of  UXO,  we  developed  a 
method  for  parametric  inversion  of  potential  UXO  targets  for  the  investigation  of  di¬ 
agnostic  parameters.  During  the  study,  we  discovered  that  data  processed  with  PCA 
can  not  be  inverted  directly  with  the  same  forward  model  used  to  invert  raw  data. 
The  forward  model  operator  must  incorporate  the  rotation  and  truncation  matrices 
used  in  the  PCA  process.  We  have  shown  that  by  incorporating  these  matrices,  the 
accuracy  of  the  recovered  diagnostic  parameters  can  be  improved. 

In  summary,  MM-1640  has  been  successful  in  understanding  the  application  of  PCA 
in  processing  TEM  data  acquired  for  UXO  applications.  We  have  developed  an  al¬ 
gorithm  ready  for  incorporation  into  standard  UXO  processing  workflows  for  the 
improvement  of  detection  of  UXO  and  the  reduction  in  false  anomalies.  We  have 
also  developed  a  consistent  parametric  inversion  algorithm  for  PCA-processed  data. 
The  ability  to  identify  and  separate  both  random  and  correlated  noise,  and  carry 
out  parametric  inversion  using  the  processed  data  is  expected  to  help  improve  the 
detection  and  discrimination  ability  and  reduce  the  cost  of  UXO  clearance. 
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